%% Script for Generating DR Fiber Coordinates
%   Written by Bryan Howell on 1/7/2013
%   Modified and commented on 2/26/2015
%   The trajectory of each DR fiber is calculated based on five reference
%   points (rostral to caudal):
%   1. Cervical (C) level 6, where modeled fibers end
%   2. most lateral DC layer where a-beta fiber trifurcates
%   3. one point that circumvents white matter
%   4. second point that circumvents white matter
%   5. a point in the foramen, where the a-beta fiber enters
%   Based on 1-5, a cubic spline is used to construct the trajectory
%
%   Anatomy
%   A-beta fibers enter the spinal cord and trifurcate. One branch projects
%   to an interneuron or a projection neuron in the spinal network encoding
%   pain (Zhang et al., 2014), one branch descends to caudal levels, and
%   one branch ascends to the dorsal column nuclei. We define a DC fiber as
%   the ascending branch of an a-beta fiber from a caudal level; and we define
%   a DR fiber as the a-beta fiber at the level of insertion plus its
%   ascending branch.

clear all;
clc;

%% Import Spinal Cord Geometry

patnum=5; % patient number (1-5)
cordgeom=load('geom_patientcord.txt'); % patient spinal cord geom.

% white matter boundary defined with an ellipse
ae=cordgeom(patnum,1); % transverse, left-right axis of ellipse
be=cordgeom(patnum,2); % transverse, down-up axis of ellipse
c_white=[0;cordgeom(patnum,3);0]; % center of ellipse

dr_th1=cordgeom(patnum,4);
dr_th2=cordgeom(patnum,5);

%% DC Fiber Geometry
%   We assume that the DC is split into 11 medial-lateral layers 
%   The lateral-most layer is where the DR fibers reside (see definitions
%   give above)

bilat = 1; % 0 = unilateral, 1 = bilateral
c_nerve = [10.5;9.5;-15]; % center of nerve at 5th ref. point

D=9; % diameter of DC fiber (ascending branch of a-beta fiber)
nodclv=10; % number of layers w/ DC fibers
axonperlv=10; % number of DC fibers per layer
noaxons=nodclv*axonperlv; % total number of fibers/axons

% I recommend savis this as a file so you can reproduce results
udist=lhsdesign(noaxons,2); % used for random sampling
%% Arc Length of Fiber
%   NoR-MYSA-FLUT-STIN-...-STIN-FLUT-MYSA-...
%   (NoR to 2nd MYSA is 1 segment)

filename=['axonarc_',num2str(D),'um.txt'];% arc length of 1 segment
axonarc=1e-3*load(filename);% um -> mm

maxarc=280; % maximum arc length
norepeats=floor(maxarc/axonarc(end)); % number of repeats (segments)

delarc=diff(axonarc)'; % distances between successive elements in 1 segment
arcval=delarc*ones(1,norepeats); % outer product 
arcval=[0;cumsum( arcval(:) )]; % string matrix from column 1 to last column
ck=(abs(arcval)<1e-10); % values less than 1x10^-10 are considered...
arcval(ck)=0; % ...zero
noarcval=length(arcval); % total number of elements

%% 1st Reference Point
%   C6 (both gracile and cuneate fasciculi present)

c6_to_t11=12; % 12 vertebral levels from C6-T11, including C6

thref_1=(64+65)/2; % halfway between 64 and 65 degress
thref_1=thref_1*pi/180; % degrees -> radians
% elliptical angle
rmag_1=ae*be./sqrt((be*cos(thref_1)).^2+(ae*sin(thref_1)).^2);
zmag_1=25*c6_to_t11; % 25 mm / level * 12 levels above insertion

% xyz coord. of 1st ref. point
xyzref_1=zeros(noaxons,3);
xyzref_1(:,1)=rmag_1*udist(:,1)*cos(thref_1);
xyzref_1(:,2)=rmag_1*udist(:,1)*sin(thref_1);
xyzref_1(:,3)=zmag_1;

%% 2nd Reference Point

% location of most lateral dermatomal level depends on patient geometry
%   - same principle as 1st ref. point (see above)
thref_2=(dr_th1+dr_th2)/2;
thref_2=thref_2*pi/180;
rmag=ae*be./sqrt((be*cos(thref_2)).^2+(ae*sin(thref_2)).^2); 
zmag=-20;

rzref_2=zeros(noaxons,2);
rzref_2(:,1)=rmag*udist(:,1);
rzref_2(:,2)=zmag*udist(:,2);

% xyz coord. of ref. point 2
xyzref_2=zeros(noaxons,3);
xyzref_2(:,1)=rzref_2(:,1)*cos(thref_2);
xyzref_2(:,2)=rzref_2(:,1)*sin(thref_2);
xyzref_2(:,3)=rzref_2(:,2);

%% 3rd Reference Point

% bording white matter (see visualization below)
thref_3=36*pi/180;
ae3=1.05*ae;
be3=1.05*be;
rmag_3=ae3*be3./sqrt((be3*cos(thref_3)).^2+(ae3*sin(thref_3)).^2);

% xyz coord. of ref. point 3
xyzref_3 = zeros(noaxons,3);
xyzref_3(:,1) = rmag_3*cos(thref_3)*ones(noaxons,1);
xyzref_3(:,2) = rmag_3*sin(thref_3)*ones(noaxons,1);
xyzref_3(:,3) = xyzref_2(:,3);

%% 4th Reference Point

% bording white matter (see visualization below)
thref_4=15*pi/180;
ae4=1.1*ae;
be4=1.1*be;
rmag_4=ae4*be4./sqrt((be4*cos(thref_4)).^2+(ae4*sin(thref_4)).^2);

% xyz coord. of ref. point 4
xyzref_4=zeros(noaxons,3);
xyzref_4(:,1)=rmag_4*cos(thref_4)*ones(noaxons,1);
xyzref_4(:,2)=rmag_4*sin(thref_4)*ones(noaxons,1);
xyzref_4(:,3)=xyzref_2(:,3);

%% 5th Reference Point

rad5=1; % diameter of nerve

% sampling points UNIFORMLY on a unit annulus
%   (if inner radius=0, unit disk)
%   Sample coordinates uniformly in polar coordinates. Transform the radii
%   so that coordinates are also uniform in Cartesian coordinates.
%   (See sampling points on a unit disk/annulus -- Wolfram Mathworld)
runit=udist(:,1);
theta=2*pi*udist(:,2);
rtrans=sqrt( runit*(rad5^2 - 0^2) + 0^2 );

% xyz coord. of ref. point 5 (all z values at 0, the level of stimulation)
xyzref_5=zeros(noaxons,3);
xyzref_5(:,2)=rtrans.*cos(theta);
xyzref_5(:,3)=rtrans.*sin(theta);

%% Translation
% translate xyz coord. to center of white matter ellipse
% - this step is only crucial if coordinating with FEM model of spine
% (see PLOS ONE article)

xyzref_1=xyzref_1+(c_white*ones(1,noaxons))';
xyzref_2=xyzref_2+(c_white*ones(1,noaxons))';
xyzref_3=xyzref_3+(c_white*ones(1,noaxons))';
xyzref_4=xyzref_4+(c_white*ones(1,noaxons))';
xyzref_5=xyzref_5+(c_nerve*ones(1,noaxons))';

%% Trajectory of DR Fiber

nodiv=1e4; % number of pieces to split the spline into
droots_coord=cell(100,1); % cell array holding coordinates

for ii=1:noaxons

    % all x, y, and z values of five reference points
    x_dr=[xyzref_5(ii,1),xyzref_4(ii,1),xyzref_3(ii,1),xyzref_2(ii,1),xyzref_1(ii,1)];
    y_dr=[xyzref_5(ii,2),xyzref_4(ii,2),xyzref_3(ii,2),xyzref_2(ii,2),xyzref_1(ii,2)];
    z_dr=[xyzref_5(ii,3),xyzref_4(ii,3),xyzref_3(ii,3),xyzref_2(ii,3),xyzref_1(ii,3)];
    xyz_dr=cat(1,x_dr,y_dr,z_dr);
    
    % ascending branch of DR fiber (will become DC of a rostral level)
    x_dc=[xyzref_2(ii,1),xyzref_1(ii,1)];
    y_dc=[xyzref_2(ii,2),xyzref_1(ii,2)];
    z_dc=[xyzref_2(ii,3),xyzref_1(ii,3)];
    xyz_dc=cat(1,x_dc,y_dc,z_dc);
    
    % fit spline to reference points
    dr_spline=cscvn(xyz_dr); % body before and after trifurcation
    dc_spline=cscvn(xyz_dc); % ascending portion only
    % ensures curvature is not too extreme
    dr_spline.coefs(10:12,:)=dc_spline.coefs(1:3,:); 
    
    so=dr_spline.breaks(1); % smallest break point (BP)
    sf=dr_spline.breaks(end); % largest BP
    sval=so:(sf-so)/nodiv:sf; % divide BP interval into small pieces 
    xyz_dr=fnval(dr_spline,sval)'; % evaluate function at given values

    % split arc into a large number of straight lines
    arc_ii=sqrt( sum(diff(xyz_dr,1,1).^2,2) );
    arc_ii=cumsum(arc_ii);

    % interpolated location - desired location
    %  - diff. > 0 = interp. value overshot
    %  - diff. < 0 = interp. value undershot
    M_arc=arc_ii*ones(1,noarcval)-( arcval*ones(1,nodiv) )';
    [junk,min_ai]=min(abs(M_arc),[],1); % closest point (row in each column)
    min_ind=sub2ind([nodiv,noarcval],min_ai,1:noarcval); % linear index in matrix
    min_a=M_arc(min_ind); % difference at closest point
    
    xyz_droot=zeros(noarcval,3);
    for jj=1:noarcval % for each desired arc length
        
        if(min_a(jj)>0)%overshot so interpolate back (rostral -> caudal)
            ra=xyz_dr(min_ai(jj)+1,:);
            rb=xyz_dr(min_ai(jj),:);
            rc=rb-ra;
        else %undershot so interpolate forward (caudal -> rostral)
            ra=xyz_dr(min_ai(jj)+1,:);
            if(min_ai(jj)<nodiv)
                rb=xyz_dr(min_ai(jj)+2,:);
                rc=rb-ra;
            else
                rb=xyz_dr(min_ai(jj),:);
                rc=ra-rb;
            end
        end
        
        % move forward/backward the difference in distance along the spline
        % (assume with enough divisions, spline in approximately linear for
        % a small division)
        xyz_droot(jj,:)=ra+rc/norm(rc)*abs( min_a(jj) );
        
    end
    
    droots_coord{ii}=xyz_droot;
    display(round(100*ii/noaxons)) % progress
    
end

% if bilateral population, assume symmetry
if(bilat==1)
   
    tmp_coord1=droots_coord;
    tmp_coord2=droots_coord;
    droots_coord=cell(2*noaxons,1);
    for k = 1:noaxons
        tmp_coord2{k}(:,1)=-tmp_coord2{k}(:,1);
        droots_coord{k}=tmp_coord1{k};
        droots_coord{noaxons + k}= tmp_coord2{k};
    end
    
    % two populations on either side
    noaxons = 2*noaxons; 
    
end

%% Storing Data
%   Condense cell array into one single matrix

m=noaxons*noarcval; % number of elements per axon * number of axons
n=3; % three spatial dimensions
Mcoord=zeros(m,n);

% condense entire cell array into a single matrix
nonodes=(noarcval-1)/8+1; % 8 elements per repeating segment
for ii =1:noaxons
    a=(ii-1)*noarcval + 1;
    b=ii*noarcval;
    Mcoord(a:b,:)=droots_coord{ii};    
end

Mcoord = Mcoord'; % transpose

file2 = ['p',num2str(patnum),'_drootxyz_',num2str(noaxons),...
         'axons_',num2str(nonodes),'nodes_',num2str(D),'um.txt'];

% % not a crucial step     
% % bounds of FEM model     
% ck1=( min( Mcoord(1,:) ) > -50 )&( max( Mcoord(1,:) ) < 50 );
% ck2=( min( Mcoord(2,:) ) > -50 )&( max( Mcoord(2,:) ) < 50 );
% ck3=( min( Mcoord(3,:) ) > -150 )&( max( Mcoord(3,:) ) < 150 );
% if(ck1&&ck2&&ck3) % if DC fiber coordinates w/in bounds of FEM model
%     dlmwrite(file2,Mcoord);
% end

%% Visualize Fibers

% unit ellipse
th_unit = 0:pi/100:2*pi;
runit   = ae*be./sqrt((be*cos(th_unit)).^2+(ae*sin(th_unit)).^2); 
xunit   = c_white(1) + runit.*cos(th_unit);
yunit   = c_white(2) + runit.*sin(th_unit);

figure;
hold on;
% reference points (minus the distal cervial reference)
plot3(xyzref_2(:,1),xyzref_2(:,2),xyzref_2(:,3),'bo');
plot3(xyzref_3(:,1),xyzref_3(:,2),xyzref_3(:,3),'bo');
plot3(xyzref_4(:,1),xyzref_4(:,2),xyzref_4(:,3),'bo');
plot3(xyzref_5(:,1),xyzref_5(:,2),xyzref_5(:,3),'bo');
% 2 ellipses to define 1 vertebral level
plot3(xunit,yunit,zeros(size(xunit)),'k-','LineWidth',2);
plot3(xunit,yunit,ones(size(xunit))-25,'k-','LineWidth',2);

for ii = 1:noaxons
    plot3(droots_coord{ii}(:,1),droots_coord{ii}(:,2),droots_coord{ii}(:,3),'k.-');

end
plot3(droots_coord{50}(:,1),droots_coord{50}(:,2),droots_coord{50}(:,3),'r.');
hold off;
title('Seeding Dorsal Roots','FontSize',36,'FontWeight','b');
xlabel('x distance (mm)','FontSize',30);
ylabel('y distance (mm)','FontSize',30);
zlabel('z distance (mm)','FontSize',30);
axis([-15,15,-5,25,-30,15]);
set(gca,'FontSize',30);
axis square;

% % XY Locations of DR Fibers
% xyloc = cat(2,xyzref_2(:,1),xyzref_2(:,2));
% xytmp = xyloc;
% xytmp(:,1) = -xytmp(:,1);
% xyloc = cat(1,xyloc,xytmp);
% plot(xyloc(:,1),xyloc(:,2),'k.');
% dlmwrite(['p',num2str(patnum),'_dr_xyloc.txt'],xyloc);
% %XYZ
% xyloc = cat(2,xyzref_2(:,1),xyzref_2(:,2),xyzref_2(:,3));
% xytmp = xyloc;
% xytmp(:,1) = -xytmp(:,1);
% xyloc = cat(1,xyloc,xytmp);
% plot3(xyloc(:,1),xyloc(:,2),xyloc(:,3),'k.');
% dlmwrite(['p',num2str(patnum),'_dr_xyzloc.txt'],xyloc);